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ABSTRACT 


Linear  anelastic  phenomena  in  wave  propagation  problems  can  be  well  modeled 
through  a  viscoelastic  mechanical  model  consisting  of  standard  linear  solids.  In  this 
paper  we  present  a  method  for  modeling  of  constant  Q  as  a  function  of  frequency 
based  on  an  explicit  closed  formula  for  calculation  of  the  parameter  fields.  The  pro¬ 
posed  method  enables  substantial  savings  in  computations  and  memory  requirements. 
Experiments  show  that  the  new  method  also  yields  higher  accuracy  in  the  modeling 
of  Q  than  e.g.,  the  Fade  approximant  method  (Day  and  Minster,  1984). 
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INTRODUCTION 


Earth  media  attenuate  and  disperse  propagating  waves.  Several  modeling  methods 
for  wave  propagation,  which  take  attenuating  and  dispersive  effects  into  account,  have 
been  presented  (e.g.,  Robertsson  et  al.  (1994),  Carcione  et  al.  (1988),  and  Martinez 
and  McMechan  (1991a))  as  have  inversion  methods  developed  for  viscoelastic  media 
(e.g.,  Blanch  and  Symes  (1994)  and  Martinez  and  McMechan  (1991b)).  Attenuating 
effects  also  have  a  large  impact  on  AVO  results  (Martinez,  1993). 

Attenuating  and  dispersive  effects  are  often  quantified  by  the  “quality  factor”,  Q. 
This  quantity  is  losely  defined  as  the  number  of  wavelengths  a  wave  can  propagate 
through  a  medium  before  its  amplitude  has  decreased  by  e~'^ .  A  more  strict  definition 
of  Q  will  be  given  in  the  following  section.  Q  has  been  found  to  be  essentially 
constant  as  a  function  of  frequency  (McDonal  et  ah,  1958)  over  the  seismic  frequency 
range,  (approximately  1-200  Hz).  Therefore  a  wave  of  higher  frequency  propagates 
a  shorter  distance  to  decay  a  certain  amount  in  amplitude,  compared  to  a  wave  of 
lower  frequency.  Futterman  (1962)  showed,  through  Kramer-Kronig’s  relation,  that  a 
constant  Q  implied  a  specific  dispersion  relation.  This  dispersion  relation  was  showed 
to  hold  for  both  shale  and  plexiglass  by  Wuenchel  (1965). 

A  correct  modeling  scheme  should  thus  yield  a  constant  Q  and  the  dispersion  re¬ 
lation  prescribed  by  Futterman  (1962).  This  is  easily  obtained  for  a  modeling  method 
based  in  the  frequency  domain,  but  is  slightly  more  cumbersome  for  modeling  meth¬ 
ods  based  in  the  time  domain.  Modeling  methods  based  in  the  frequency  domain 
are  essentially  limited  to  1-D  or  layered  (through  propagator  matrices)  media  and 
we  will  therefore  focus  on  “time  domain  methods” .  An  approximately  constant  Q  is 
usually  constructed,  by  using  simple  viscoelastic  building  blocks,  to  obtain  a  numer¬ 
ically  efficient  modeling  method  based  in  the  time  domain.  Day  and  Minster  (1984) 
showed  that  the  best  Fade  approximant  to  a  constant  Q  is  a  so-called  Standard  Lin¬ 
ear  Solid  (SLS).  Day  and  Minster  (1984)  also  described  a  method  to  connect  several 
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SLSs  in  parallel  to  yield  an  excellent  approximation  to  a  constant  Q  in  a  pre-defined 
frequency  band.  The  method  requires  several  SLSs  to  yield  a  good  result  and  each 
SLS  is  described  by  an  individual  parameter  set.  These  two  drawbacks  makes  the 
approximation  method  expensive,  with  respect  to  both  memory  and  time,  to  use  for 
viscoelastic  wavepropagation  modeling  methods.  The  method  is  also  quite  compli¬ 
cated,  since  it  for  instance  is  necessary  to  find  roots  of  Legendre  polynomials. 

In  this  paper  we  present  a  method  to  find  approximations  to  a  constant  Q  which  is 
easier  to  use  than  Day  and  Minster’s  (1984)  method  and  yields  good  approximations 
using  only  a  few  SLSs.  The  viscoelastic  description,  used  in  the  method,  yields  savings 
during  the  viscoelastic  wave  simulation  as  well.  We  have  chosen  to  call  the  method, 
the  “r-method”,  after  the  notation  used  in  deriving  the  algorithm. 

We  will  start  by  reviewing  some  basic  viscoelastic  concepts  and  definitions.  We 
will  then  continue  by  describing  our  method  and  finally  compare  results  from  finite- 
difference  viscoelastic  wave  propagation  simulations,  where  the  attenuation  model 
was  found  using  the  r-method,  a  single  SLS  approximation,  the  “Fade  approximant” 
method  (Day  and  Minster,  1984),  and  an  essentially  analytic  solution  based  on  results 
from  McDonal  et  al.  (1958),  Futterman  (1962),  and  Wuenschel  (1965).  The  analytic 
solution  is  obtained  through  a  frequency  domain  method  similar  to  the  method  de¬ 
scribed  by  Martinez  and  McMechan  (1991a). 

LINEAR  VISCOELASTIC  MODEL 

The  constitutive  relation  for  a  linear  elastic  2-D  {i,j,k  =  x,y)  or  3-D  — 

x,y,z)  homogeneous  solid  is, 

^ij  —  T  (1) 

where  cr^j  are  the  components  of  the  stress  tensor,  A  and  y  the  Lame  constants, 
and  Cij  the  components  of  the  strain  tensor.  For  a  linear  viscoelastic  material  the 
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multiplications  in  equation  (1)  become  convolutions  and  the  Lame  constants  time 
dependent,  so  that, 


A  *  EkkSij  +  2M  *  £ 


(2) 


(Christensen,  1982).  The  dot  denotes  derivative  in  time.  Thus,  a  current  deformation 
state  is  related  to  previous  ones  through  the  stress  relaxation  functions,  A(t)  and  M(t). 
It  is  possible  to  express  the  time  derivative  of  Sij  as. 


£ij 


(3) 


where  Vi  is  the  velocity.  Let  us  define, 


n  =  A  +  2M.  (4) 

Equations  (2),  (3),  and  (4)  then  yields, 

&ij  =  (il  —  2M)  ♦  dkVk  +  M  *  diVj.  (5) 

Together  with  Newton’s  second  law, 

Vi  =  -djaij,  (6) 

P 

equation  (5)  describes  viscoelastic  wave  propagation.  The  definition  of  n(f)  allows 
us  to  independently  define  quality  factors,  Q,  for  P-  and  S'-waves.  Robertsson  et  al. 
(1994)  describe  an  efficient  0(2,4)  finite-difference  scheme  to  solve  the  problem. 

The  SLS  (a  spring  and  a  dashpot  in  series  in  parallel  with  a  spring)  has  been 
shown  to  be  a  quite  general  mechanical  viscoelastic  model  (e.g..  Day  and  Minster 
(1984),  Pipkin  (1986),  and  Blanch  et  al.  (1993)).  An  array  of  L  SLS  has  the  stress 
relaxation  function, 

r(i)  =  (^1  -  2  (i  -  :^)  (7) 

where  0{t)  is  the  Heaviside  function,  Mr  is  the  relaxed  stress  modulus  corresponding 
to  r(f)  (Pipkin,  1986),  and  and  are  the  stress  and  strain  relaxation  times  for 
the  /th  SLS  (Pipkin,  1986).  With  our  formulation  we  therefore  have. 
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M(()  =  f  -  E  (l  -  “)  "■""")  «(0.  (9) 

where  tt  =  A  +  2/r,  and  rh  and  are  the  strain  relaxation  times  for  P-  and  5- waves 
for  the  /th  SLS.  The  same  stress  relaxation  times,  r^/,  may  be  used  for  both  P-  and 
S'-waves.  In  this  paper,  we  will  also  refer  to  an  SLS  as  a  relaxation  mechanism. 

The  complex  stress  modulus,  Mc{eo),  is  defined  as  the  Fourier  transform  of  the 
stress  relaxation  function  (Pipkin,  1986).  The  quality  factor,  Q,  is  defined  as, 

(Bourbie  et  ah,  1987).  Equation  (10)  defines  Q  as  the  number  of  wavelengths  a  pulse 
may  propagate  before  its  amplitude  drops  by  a  factor  of  e~'"  (White,  1992).  Q  is  thus 
generally  a  function  of  frequency.  For  an  array  of  standard  linear  solids  equations  (7) 
and  (10)  yields. 


^  Uj{Tel  -  T„l) 


We  may  thus  express  the  viscoelastic  medium  in  terms  of  two  arrays  of  SLS,  n(f) 
and  M(t)  (equations  [8]  and  [9]),  given  the  Q  vs.  frequency  relations  for  P-  and 
S'-waves,  the  P-  and  S-velocities,  and  the  density  of  the  medium.  Since  attenuation 
and  dispersion  are  related  through  Kramer-Kronig’s  relation  (Futterman,  1962),  the 
choice  of  a  physically  realistic  Q  vs.  frequency  relation  yields  an  equally  physical 
dispersion  relation.  We  will  now  derive  a  formula  on  closed  form  to  calculate  the  t^i 
and  the  in  equation  (11)  for  a  desired  constant  Q,  for  P-  and  S-waves  respectively. 

THE  r-METHOD 


The  basic  arguments  for  the  r-method  are. 
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•  simplicity  in  determining  a  constant  Q  model, 

•  memory  savings  and  less  calculations  in  forward  modeling  of  viscoelastic  waves, 

•  the  convenience  in  having  only  one  variable  describing  the  magnitude  of  disper¬ 
sion  and  attenuation  in  inversion  routines. 

We  will  later  see  that  the  r-method  also  yields  more  accurate  results  at  a  low  compu¬ 
tational  cost  than  e.g.,  the  Fade  approximant  method.  We  will  assume  that  we  want 
to  approximate  a  constant  Q  in  a  pre-defined  frequency  interval  and  that  there  exist 
L  preset  stress  relaxation  times,  (see  below). 

The  r-method  is  based  on  the  simple  observation  that  the  level  of  attenuation 
caused  by  a  SLS  can  be  determined  by  a  dimensionless  (frequency  scale-independent) 
variable,  r.  If  we  define  r  as, 

^  _  '^el  '^crl 

'^crl  Tt/ 

the  inverse  of  Q  for  one  SLS  can  be  written  as, 

1  +u}'^T^{l  -fr)’ 

It  is  easily  seen  that  the  behavior  in  frequency  is  essentially  determined  by  To-  and 
the  magnitude  of  Q  for  the  SLS  essentially  by  r.  This  is  not  true  for  r  larger  than,  or 
even  close  to,  1.  As  we  shall  see,  r  <C  1  generally.  We  have  plotted  Q  as  a  function  of 
frequency  for  a  SLS  with  two  diffferent  and  two  different  r  in  Figure  1,  to  illustrate 
the  effect  of  the  different  parameters.  The  magnitude  is  the  same  but  the  curve  is 
displaced  in  frequency  when  is  changed  but  r  remains  constant.  The  magnitude 
of  the  curve  changes  with  r  on  the  other  hand,  but  the  curve  does  not  move  in 
frequency  when  is  held  constant.  We  may  therefore  use  r  as  a  means  to  determine 
the  magnitude  for  several  different  SLSs,  which  assume  maximum  attenuation  at 
different  frequencies.  For  one  SLS  r  where  Qmin  is  the  global  minimum  for 

Q  as  a  function  of  frequency.  For  more  than  one  SLS,  r  is  even  smaller.  Q  is  rarely 
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less  than  20  for  any  real  materials,  implying  that  r  <C  1.  Using  the  parameter  t  to 
tune  an  array  of  SLSs,  and  assuming  that  1  +  r  srf  1,  equation  (11)  yields, 


Ui 

/=!  ^ 


+  u;2r,V 


(14) 


In  this  (approximate)  expression  Q~^  is  linear  in  r.  We  can  therefore  easily  find  the 
best  approximation  in  the  least  squares  sense  over  a  pre-defined  frequency  range  to 
any  Qo,  by  minimizing  over  r  the  expression. 


(Q  \^,T<Thr)  -  duj.  (15) 

Juja 

To  find  the  minimum  we  set  the  derivative  of  J  with  respect  to  r  to  zero,  and  solve 
for  T. 


-  =  (Q  Kr.,.r)  -  Q„  ) - - - du,  =  0 


(16) 


If  we  define, 


E 


1=1 


we  can  write  equation  (16)  as. 


(17) 


dj  o 

~  =  2  T  {F{UJ,  T„i))  -  Qq  ^F{uj,  T^i)  du  =  0. 

UT  J  UJa 

We  can  now  solve  for  r,  and  find. 


(18) 


r  = 


'  F{u,T„i)du; 
_ _ 

^0  /  {F{u>,T„i)f  du 


(19) 


The  two  integrals  in  equation  (19)  can  be  solved  analytically,  which  yields  an 
explicit  formula  for  r.  The  integrands  (or  methods)  can  be  found  in  any  standard 
work  on  integral  calculus  or  collection  of  mathematical  formulas  (e.g.,  Rade  and 
Westergren  (1988)).  It  is  of  course  possible  to  solve  the  integrals  by  some  numerical 
method  as  well.  The  final  formula  for  r  is. 
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T  = 


1=1 


(20) 


where 


Qo 


^  7i/  +  2  ^  ^  hki 


1=1  1=1  k=i 
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Ioi  = 

hi  = 
hik  = 


1 


1 

'^al'^crk 

-2  _  ,2 
^ak  ^erl 


2  2 


)1 


log  (l  + 
arctan  (uiTai)  — 


W6 

coa 


^Tal 


‘^b 


1+U,2r2_  _ 

arctan  {u>Tt,i )  arctan  (coTf^k) 


Tel 


'^ak 


(21) 

(22) 

(23) 


The  r-method  generally  yields  a  slightly  too  small  r,  i.e.,  slightly  too  large  Q, 
caused  by  the  approximation  r  -C  1.  By  using  a  slightly  smaller  Qq  in  equation  (20) 
than  desired,  it  is  possible  to  find  the  r  which  yields  the  desired  Q.  The  effect  of 
underestimating  r  is  more  pronounced  for  a  larger  number  of  relaxation  mechanisms 
and  a  lower  Q.  An  example  of  the  r-method  is  shown  in  Figure  3,  where  five  relaxation 
mechanisms  have  been  used  to  find  a  Q  approximately  equal  to  20.  Hence  we  used  a 
Qq  approximately  equal  to  18  in  equation  (20)  to  calculate  r. 

If  several  r  are  to  be  calculated  for  an  identical  distribution  and  frequency 
intervals,  but  different  Qo,  it  suffices  to  calculate  the  integrals  /q/,  /k,  and  l2ik  and 
corresponding  the  sums  once.  This  can  be  seen  from  equation  (20)  since  the  integrals 
are  independent  of  Qq.  To  determine  a  new  r  it  is  thus  sufficient  to  divide  by  a  new 
Qo- 

The  distribution  of  the  t^i  is  not  determined  by  the  r-method.  Distributing  the 
Tel  logarithmically  over  the  frequency  range  of  interest  generally  yields  good  approxi¬ 
mations  to  a  constant  Q.  This  is  r^/  =  1/cu/,  about  one  per  one-two  octaves,  which  is 
a  good  ad  hoc  rule  of  thumb  (cf.,  Liu  et  al.  (1976),  Blanch  et  al.  (1993),  and  example 
in  Figure  3  and  Table  4).  A  closed  formula  for  the  optimal  distribution  of  the  Tei  is 
still  an  open  question. 


The  use  of  a  single  parameter  determining  the  magnitude  of  several  SLSs  also 
yields  considerable  savings  for  forward  modeling  of  viscoelastic  waves.  There  are 
both  less  memory  and  calculations  required  for  a  certain  problem  (see  section  Imple¬ 
mentation  Considerations). 

NUMERICAL  EXPERIMENTS 

In  this  section  we  benchmark  the  r-method  against  the  Fade  approximant  method 
(Day  and  Minster,  1984),  a  single  relaxation  mechanism  (one  SLS)  approximation, 
and  an  analytical  solution  which  we  will  start  by  deriving.  We  compare  seismograms 
from  1-D  finite-difference  simulations  using  the  method  described  by  Robertsson  et 
al.  (1994).  The  results  can  immediately  be  generalized  to  2-D  and  3-D,  since  P-  and 
S'-waves  are  modeled  independently.  We  use  models  with  different  attenuation  as  well 
as  (^-approximations  based  on  different  numbers  of  relaxation  mechansims  (or  SLS) 
to  evaluate  the  different  methods. 

An  analytical  solution  for  anelastic  wave  propagation  in  a  homogeneous  medium 
may  be  obtained  through  the  principle  of  correspondence  (Bland,  1960).  From 
D’Alembert’s  solution  we  obtain  the  impulse  response  for  the  one-dimenstional  case 
as, 

<^(u;,x)  =  (24) 

where  i  —  a/— 1,  x  is  the  offset  from  the  initial  location  of  the  source,  and, 

n(u;)  =  c{u)  ^1  -f 

where  sgn{u)  =  —1  for  a?  <  0,  and  sgn{Lo)  =  1  for  a;  >  0.  c{uj)  and  Q{io)  are  the 
desired  velocity  and  Q  relations.  An  arbitrary  source  wavelet  can  thus  be  convolved 
with  the  impulse  response  in  equation  (24)  to  yield  the  solution  as  a  function  of  offset 
and  time,  enabling  the  calculation  of  analytic  seismograms  for  desired  velocity  and 
(^-relations  as  functions  of  frequency. 
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Futterman  (1962)  proposed  three  different  analytical  expressions  to  model  as  a 
function  of  frequency  in  real  earth  media  based  on  theoretical  rock  mechanics  studies. 
The  velocity  vs.  frequency  relation  is  given  from  the  Q  vs.  frequency  relation  through 
Kramer-Kronig’s  relation.  Martinez  and  McMechan  (1991a)  chose  one  of  the  expres¬ 
sions  to  obtain  a  r-p  domain  solution  for  anelastic  wave  propagation.  Using  shale  and 
Plexiglass,  Wuenschel  (1965)  showed  that  another  of  Futtermans  (1962)  expressions 
fits  very  close  to  experimental  results.  This  is  the  one  we  used  to  benchmark  the 
different  (^-modeling  algorithms.  The  velocity  is  given  by, 

c{uj)  =  Co/  ^1  -  I  (t^/^o)^  ~  ^ 

and  thus  Q  is, 

Q{u)  =  Qo  I  (W^o)^  ”  ^ 

tuo  is  a  reference  frequency,  much  smaller  than  the  frequencies  modeled.  In  our  ex¬ 
periments  we  used  a;o=0.01  rad/s,  but  the  value  of  uq  is  not  crucial  for  the  model,  cq 
is  the  reference  velocity  chosen  to  yield  the  desired  velocity  at  the  center  frequency  of 
the  wavelet  in  the  simulations,  and  Qq  is  the  analogous  reference  Q-value.  In  our  sim¬ 
ulations  we  will  refer  to  the  velocity  and  the  Q  as  the  values  at  the  center  frequency 
of  the  wavelet  used  in  the  experiments  (10  Hz). 

Test  case 

The  source  in  the  simulations  is  a  Ricker  wavelet  with  an  amplitude  of  1  and  a 
center  frequency  of  10  Hz,  which  yields  a  band  width  between  2  and  25  Hz  for  a  40  dB 
threshhold.  Number  of  gridpoints  and  timestep  for  the  viscoelastic  simulation  were 
chosen  in  accordance  with  rules  outlined  by  Robertsson  et  al.  (1994).  The  density, 
p,  is  1000  kg/m^,  and  Co=1500  m/s.  We  chose  two  models.  In  the  first  we  used  a 
moderate  value  of  Qo=100.  This  is  typical  for  shallow  rocks  or  consolidated  sediments 
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(e.g.,  Bourbie  et  al.  (1987)  and  Hamilton  (1980)).  In  the  second  model  Qo=20,  which 
is  a  low  value  characteristic  of  shallow  unconsolidated  sediments  on  the  seafloor  (e.g., 
Hamilton  (1980)).  The  receivers  are  located  at  3000  m  and  6000  m  from  the  source, 
corresponding  to  20  and  40  wavelengths  for  the  center  frequency  of  the  wavelet. 

Q-modeling. — We  investigated  two  cases.  First,  we  employed  five  relaxation 
mechanisms  to  approximate  the  Q=20  model.  For  strong  attenuation,  accurate  mod¬ 
eling  of  Q  is  naturally  more  crucial  than  for  cases  where  the  attenuation  is  weaker. 
We  do  not  show  results  from  approximations  of  a  Q  of  100  using  five  relaxation 
mechanisms  since  the  excellent  results  we  achived  for  the  Q=20  approximation  al¬ 
ready  illustrates  the  case  of  five  relaxation  mechanisms.  Second,  we  employed  two 
relaxation  mechanisms  to  approximate  both  the  Q=20  and  the  <5=100  models.  As 
described  above,  we  have  constrained  the  distribution  of  so  that  for  each  of  the 
Q-optimization  methods,  the  Tai  are  identical  both  for  <5=100  and  <5=20.  This  yields 
optimally  low  memory  requirements,  which  is  particularly  important  in  2-D  and  3-D 
applications.  The  relaxation  times  obtained  from  the  different  <5-optimization  meth¬ 
ods  are  listed  in  Tables  1,  2,  3,  4,  and  5. 

Equation  (11)  may  be  used  to  calculate  <5  as  a  function  of  frequency  for  the  dif¬ 
ferent  <5-optimization  methods.  In  Figure  2  we  have  plotted  these  curves  when  using 
two  relaxation  mechanisms  to  approximate  a  <5  of  20  between  2  and  25  Hz.  As  we  can 
see  the  r-method  yields  a  close  to  constant  Q  over  the  entire  frequency  band.  The 
single  relaxation  mechanism  approximation  yields  a  Q  that  roughly  varies  between 
18  and  30.  The  chosen  source  wavelet  has  its  energy  concentrated  around  10  Hz,  and 
since  the  single  relaxation  mechanism  approximation  varies  between  a  <5  of  18  and 
22  between  5  and  20  Hz,  we  can  expect  fairly  good  simulation  results  also  for  this 
optimally  computational  inexpensive  <5-modeling  method.  The  Fade  approximant 
method  yields  a  quite  poor  result,  barely  better  than  the  single  relaxation  mechanism 
approximation.  The  Fade  approximant  method  is  in  some  sense  optimal,  but  con- 
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verges  slowly  as  the  number  of  approximants  are  increased,  when  a  constant- Q  over 
a  relatively  short  frequency  band  is  desired. 

In  Figure  3  we  show  the  (^-approximations  when  using  five  mechanisms  to  approx¬ 
imate  a  (5  of  20  between  2  and  25  Hz.  The  r-method  yields  very  close  to  constant 
Q  over  the  frequency  band.  Also,  the  Fade  approximant  method  now  gives  a  much 
better  approximation  compared  to  the  case  when  using  only  two  approximants.  It  is 
however  still  worse  than  the  approximation  obtained  through  the  r-method. 

In  Figure  4  we  illustrate  the  (^-approximations  when  using  two  mechanisms  to 
approximate  a.  Q  oi  100  between  2  and  25  Hz.  Again,  the  r-method  yields  the  best 
result  and  the  Fade  approximant  method  yields  a  result  similar  to  that  of  the  single 
relaxation  mechanism  approximation. 

Simulations. — In  Figure  5  we  have  plotted  the  seismograms  collected  at  6000  m 
offset  when  approximating  a  constant  Q  of  100  with  two  relaxation  mechanisms.  As 
we  can  see  all  methods  reproduce  the  analytic  solution  fairly  well,  although  the  r- 
method  clearly  is  the  most  successful.  The  Fade  approximant  method  and  the  single 
relaxation  mechanism  approximation  both  yield  worse  but  reasonably  accurate  and 
similar  solutions.  The  latter  is  not  a  surprise  since  the  corresponding  Q  vs.  frequency 
relations  are  very  similar  (see  Figure  4). 

In  Figure  6  we  show  the  seismograms  collected  at  3000  m  offset  when  approxi¬ 
mating  a  constant  Q  of  20  using  two  relaxation  mechanisms.  All  methods  display 
highly  accurate  results,  even  though  the  r-method  again  is  most  successful.  In  Fig¬ 
ure  7  we  show  the  analogous  seismograms  collected  at  6000  m  offset.  Notice  that 
the  amplitude  of  the  analytic  solution  now  only  is  approximately  1  percent  of  the 
initial  amplitude.  Furthermore,  for  40  wavelengths  path  of  propagation,  cylindrical 
and  spherical  spreading  would  significantly  contribute  to  even  stronger  decay  of  the 
waveform  in  2-D  and  3-D  simulations.  There  are  few  seismic  applications,  if  any, 
where  accurate  modeling  of  such  extreme  cases  is  of  interest.  We  achive  a  highly 
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accurate  solution  using  the  r-method  with  only  two  relaxation  mechanisms.  The 
Fade  approximant  method  as  well  as  the  single  relaxation  mechanism  approximation 
yield  less  accurate  results.  The  single  relaxation  mechanism  is  however  the  trivial 
and  always  least  computationally  expensive  constant  Q  approximation.  We  find  it 
remarkable  that  such  a  simple  approximation  yields  a  for  most  applications  suffi¬ 
ciently  accurate  result  after  40  wavelengths  propagation  of  a  Ricker  wavelet  through 
a  medium  with  a  Q  of  20.  Since  most  high  frequency  components  are  attenuated, 
the  seismograms  we  collect  at  6000  m  offset  is  reflected  by  the  character  of  the  Q 
vs.  frequency  curves  at  lower  frequencies  (somewhere  around  5  Hz).  These  compo¬ 
nents  were  much  weaker  in  the  original  wavelet  but  are  not  subjected  to  as  strong 
attenuation  due  to  the  fewer  wavelengths  in  the  path  of  propagation. 

Finally  we  show  the  results  when  using  five  relaxation  mechanisms  to  approxi¬ 
mate  a  constant  Q  of  20.  At  3000  m  offset  (Figure  8),  the  analytic  solution  is  well 
reproduced  by  all  methods,  as  was  the  case  in  Figure  6  when  using  only  two  relax¬ 
ation  mechanisms.  Again,  the  r-method  yields  the  most  accurate  reult.  In  Figure  9 
we  show  seismograms  collected  at  6000  m  offset.  The  r-method  yields  an  excellent 
result;  The  collected  seismogram  is  very  close  to  the  analytic  solution  for  this  the 
most  extreme  case  of  strong  attenuation  and  long  path  of  propagation.  As  is  evident 
from  Figure  3,  the  method  models  wave  propagation  just  as  well  over  the  entire  fre¬ 
quency  band  between  2  and  25  Hz.  The  practically  inconvenient  and  expensive  Fade 
approximant  method  yields  a  quite  poor  result,  only  marginally  more  accurate  than 
the  single  relaxation  mechanism  approximation. 

IMPLEMENTATION  CONSIDERATIONS 

Modelling  of  viscoelastic  wave  propagation  has  been  regarded  as  considerably  more 
computationally  expensive  than  purely  elcistic  wave  propagation.  Robertsson  et  al. 
(1994)  presented  a  highly  efficent  viscoelastic  finite- difference  scheme  enabling  mod- 
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eling  of  widely  varying  Q  in  highly  heterogeneous  media.  2-D  as  well  as  3-D  schemes 
have  been  implemented  at  a  computational  cost  that  only  marginally  exceeds  that  of 
analogous  purely  elastic  finite-difference  schemes  (Robertsson  et  ah,  1994).  Compu¬ 
tational  efficiency  may  be  stated  in  terms  of  the  number  of  parameter  and  variable 
fields  that  are  required  to  be  stored  simultaneously  and  the  number  of  calculations 
per  grid-point  and  time  step  required  for  a  particular  scheme.  In  Table  6  we  compare 
the  computational  efficiency  for  2-D  and  3-D  viscoelastic  schemes  using  the  different 
Q-modeling  algorithms  to  the  analogous  elastic  schemes.  Both  the  memory  require¬ 
ments  and  the  number  of  calculations  are  very  much  dependent  on  the  Q-modeling 
algorithm  employed.  From  Table  6  we  conclude  that  the  T-method  is  superior  to 
the  Fade  approximant  method  also  when  it  comes  to  computational  efficiency.  The 
simplest  Q-modeling  method,  the  single  relaxation  mechanism  approximation,  is  only 
marginally  more  expensive  than  purely  elastic  modelling. 

As  we  have  seen  the  single  relaxation  mechanism  approximation  yields  results  that 
are  sufficiently  accurate  for  most  practical  purposes.  Q  is  a  physical  quantity  that  is 
difficult  to  measure  with  high  accuracy  in  real  earth  media  (White,  1992).  Also,  it 
is  still  debated  how  realistic  a  constant  as  a  function  of  frequency  model  always 
is  (Bourbie  et  ah,  1987).  Thus,  the  uncertainty  of  the  modeling  may  very  well  be 
greater  than  the  inaccuracy  of  a  single  relaxation  mechansim  constant  Q  approxima¬ 
tion.  There  is  however  another  reason  for  using  more  than  one  relaxation  mechanisms. 
Absorbing  boundaries  can  be  implemented  by  decreasing  the  attenuation  to  as  low 
as  Q=2  within  a  narrow  frame  around  the  desired  grid  (Robertsson  et  ah,  1994).  To 
avoid  impedance  reflections  it  is  important  to  tune  the  Lame  parameters  for  constant 
velocity.  The  deviation  from  a  constant  Q  for  such  low  values  of  Q  is  greater  than  for 
the  more  physically  realistic  values  of  Q  within  the  desired  grid,  leading  to  velocity 
variations  within  the  frequency  band  of  interest  in  the  absorbing  frame.  Compensa¬ 
tion  may  thus  only  be  achived  by  employing  additional  relaxation  mechanisms  for  a 
closer  to  constant  Q  in  the  absorbing  region.  For  a  wavelet  with  a  frequency  content 
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similar  to  that  of  the  Ricker  wavelet,  we  found  that  two  relaxation  mechanisms  yield 
sufficiently  constant  Q  in  the  absorbing  frame. 

CONCLUSIONS 

We  have  outlined  the  methodology  of  an  algorithm,  the  r-method,  to  model  con¬ 
stant  Q  vs.  frequency  for  viscoelastic/ viscoacoustic  wave  propagation  simulators.  The 
T-method  enables  considerable  savings  of  memory  and  computations  for  both  2-D  and 
3-D  viscoelastic  finite-difference  schemes  and  is  much  easier  to  implement  compared 
to  other  methods,  such  as  the  Fade  apprOximant  method  (Day  and  Minster,  1984) 

We  performed  a  series  of  1-D  experiments  to  compare  the  r-method  to  an  analyti¬ 
cal  solution  based  on  the  principle  of  correspondence  (Bland,  1960)  and  a  constant-Q 
model  by  Futterman  (1962).  We  also  compared  the  method  to  the  Fade  approxi- 
mant  method  and  a  single  relaxation  mechanism  approximation.  We  found  that  the 
r-method  also  yields  the  highest  accuracy  in  the  simulations.  The  Fade  approximant 
method  should  however  yield  the  most  accurate  results  for  cases  with  large  numbers 
(probably  at  least  10)  of  SLSs  (Day  and  Minster,  1984). 

The  single  relaxation  mechanism  approximation  produced  reasonably  accurate 
results.  For  most  applications  this  method  yields  sufficiently  accurate  modeling  of  Q. 
In  practice,  we  see  no  reason  for  employing  more  accurate  approximations  than  the  r- 
method  using  two  relaxation  mechanisms.  This  method  also  enables  implementation 
of  highly  efficient  absorbing  boundaries. 
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TABLES 


1 

Tal  (ms) 

I 

99.472 

I 

7.2343 

TABLE  1.  Optimized  for  two  relaxation  mechanisms  to  yield  a  constant  Q  between 


2  and  25  Hz,  using  the  r-method.  r=0. 10110  for  <5=20,  and  r=0. 019156  for  Q=100. 


I 

Tal  (rns) 

Tei  (ms)  ((9=20) 

T^i  (ms)  ((9=100) 

B 

5.7653 

5.9428 

5.7946 

1 

21.495 

23.471 

21.882 

TABLE  2.  Optimized  relaxation  mechanisms  (for  a  total  of  two)  to  yield  a  constant  Q 
between  2  and  25  Hz,  using  the  Fade  approximant  method. 
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(ms) 

Ts  (ms)  ((5=20) 

Te  (ms)  ((5  =  100) 

15.056 

16.824 

15.374 

TABLE  3.  Optimized  relaxation  mechanism  to  yield  a  constant  Q  between  2  and  25  Hz, 
for  a  single  relaxation  mechanism. 


I 

T<rl  (ms) 

fl 

265.26 

B 

52.203 

B 

10.273 

B 

2.0218 

B 

0.39789 

TABLE  4.  Optimized  for  five  relaxation  mechanisms  to  yield  a  constant  Q  between 
2  and  25  Hz,  using  the  r-method.  r=:0. 060168  for  Q=20. 
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T„i  (ms) 

Tsi  (ms)  ((3=20) 

1 

1.1132 

1.1265 

B 

1.3792 

1.4033 

B 

2.1214 

2.1748 

B 

4.5928 

4.7730 

B 

22.466 

24.341 

TABLE  5.  Optimized  relaxation  mechanisms  (for  a  total  of  five)  to  yield  a  constant  Q 
between  2  and  25  Hz,  using  the  Fade  approximant  method. 
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Method 

Memory  (2-D) 

Calculations  (2-D) 

Memory  (3-D) 

Calculations  (3-D) 

Single  mech. 

13 

120 

20 

220 

T  2  mech. 

16 

130 

26 

240 

r  5  mech. 

25 

170 

44 

310 

Fade  2  mech. 

18 

140 

28 

255 

Fade  5  mech. 

33 

200 

52 

360 

Elastic 

8 

70 

12 

160 

TABLE  6.  The  computational  efficiency  of  viscoelastic  finite-difference  schemes  as 
(Robertsson  et  ah,  1994)  implemented  them  when  employing  the  different  Q-modeUng  meth¬ 
ods,  compared  to  the  analogous  purely  elastic  finite-difference  schemes.  The  computational 
efficiency  is  characterized  by  the  number  of  parameter  or  variable  fields  required  to  be  stored 
simultaneously  and  by  the  number  of  calculation  per  gridpoint  and  time  step  necessary. 
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FIGURES 

FIG.  1.  Q  as  a  function  of  frequency  for  two  different  r  and  two  different 
Solid:  T  =  4.6212  x  10“^,  =  1.5915  x  10“^  corresponds  to  10  Hz.  Dashed: 

r  =  4.6212  X  10“^,  =  1.5915  x  10“^  corresponds  to  100  Hz.  Dash-dotted: 

r  =  2.3106  X  10“^,  T„  —  1.5915  x  10“^  corresponds  to  10  Hz. 

FIG.  2.  Approximations  to  a  constant  Q  of  20  between  2  and  25  Hz.  Qq  18 
in  algorithm.  Solid:  Desired  Q.  Dashed:  r-method  using  two  relaxation  mechanisms. 
Dotted:  Fade  approximant  method  using  two  relaxation  mechanisms.  Dash-dotted: 
Single  relaxation  mechanism  approximant. 

FIG.  3.  Approximations  to  a  constant  Q  of  20  between  2  and  25  Hz.  Qq  ^  18 
in  algorithm.  Solid:  Desired  Q.  Dashed:  r-method  using  five  relaxation  mechanisms. 
Dotted:  Fade  approximant  method  using  five  relaxation  mechanisms.  Dash-dotted: 
Single  relaxation  mechanism  approximant. 

FIG.  4.  Approximations  to  a  constant  Q  of  100  between  2  and  25  Hz.  Qq  ~  95 
in  algorithm.  Solid:  Desired  Q.  Dashed:  r-method  using  two  relaxation  mechanisms. 
Dotted:  Fade  approximant  method  using  two  relaxation  mechanisms.  Dash-dotted: 
Single  relaxation  mechanism  approximant. 

FIG.  5.  Seismograms  collected  at  6000  m  offset  in  the  model  with  a  constant  Q 
of  approximation  of  100  between  2  and  25  Hz.  Solid:  Analytical  solution.  Dashed:  r- 
method  using  two  relaxation  mechanisms.  Dotted:  Fade  approximant  method  using 
two  relaxation  mechanisms.  Dash-dotted:  Single  relaxation  mechanism  approximant. 

FIG.  6.  Seismograms  collected  at  3000  m  offset  in  the  model  with  a  constant  Q 
of  approximation  of  20  between  2  and  25  Hz.  Solid:  Analytical  solution.  Dashed:  r- 
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method  using  two  relaxation  mechanisms.  Dotted:  Fade  approximant  method  using 
two  relaxation  mechanisms.  Dash-dotted:  Single  relaxation  mechanism  approximant. 

FIG.  7.  Seismograms  collected  at  6000  m  offset  in  the  model  with  a  constant  Q 
of  approximation  of  20  between  2  and  25  Hz.  Solid:  Analytical  solution.  Dashed:  r- 
method  using  two  relaxation  mechanisms.  Dotted:  Fade  approximant  method  using 
two  relaxation  mechanisms.  Dash-dotted:  Single  relaxation  mechanism  approximant. 

FIG.  8.  Seismograms  collected  at  3000  m  offset  in  the  model  with  a  constant  Q 
of  approximation  of  20  between  2  and  25  Hz.  Solid:  Analytical  solution.  Dashed:  r- 
method  using  five  relaxation  mechanisms.  Dotted:  Fade  approximant  method  using 
five  relaxation  mechanisms.  Dash-dotted:  Single  relaxation  mechanism  approximant. 

FIG.  9.  Seismograms  collected  at  6000  m  offset  in  the  model  with  a  constant  Q 
of  approximation  of  20  between  2  and  25  Hz.  Solid:  Analytical  solution.  Dashed:  r- 
method  using  five  relaxation  mechanisms.  Dotted:  Fade  approximant  method  using 
five  relaxation  mechanisms.  Dash-dotted:  Single  relaxation  mechanism  approximant. 
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FIG.  1.  Q  as  a  function  of  frequency  for  two  different  r  and  two  different  Solid: 
r  =  4.6212  X  10“^,  =  1.5915  X  10“^  corresponds  to  10  Hz.  Dashed:  r  =  4.6212  X  lO”^, 

T„  =  1.5915  X  10“^  corresponds  to  100  Hz.  Dash-dotted:  t  =  2.3106  x  10“^, 
Ter  =  1.5915  X  10“^  corresponds  to  10  Hz. 
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f  (Hz) 

FIG.  2.  Approximations  to  a  constant  Q  of  20  between  2  and  25  Hz.  Qo  ~  18  in 


algorithm.  Solid:  Desired  Q.  Dashed:  r-method  using  two  relaxation  mechanisms.  Dot¬ 
ted:  Fade  approximant  method  using  two  relaxation  mechanisms.  Dash-dotted:  Single 
relaxation  mechanism  approximant. 
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FIG.  3.  Approximations  to  a  constant  Q  of  20  between  2  and  25  Hz.  Qo  w  18  in  al 
gorithm.  Solid:  Desired  Q.  Dashed:  r-method  using  five  relaxation  mechanisms.  Dotted 
Fade  approximant  method  using  five  relaxation  mechanisms.  Dash-dotted:  Single  relaix 
ation  mechanism  approximant. 
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FIG.  4.  Approximations  to  a  constant  Q  of  100  between  2  and  25  Hz.  Qo  ~  95  in 
algorithm.  Solid:  Desired  Q.  Dashed:  r-method  using  two  relaxation  mechanisms.  Dot¬ 
ted:  Fade  approximant  method  using  two  relaxation  mechanisms.  Dash-dotted:  Single 
relaxation  mechanism  approximant. 
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t(s) 

FIG.  5.  Seismograms  collected  at  6000  m  offset  in  the  model  with  a  constant  Q  of 
approximation  of  100  between  2  and  25  Hz.  Solid:  Analytical  solution.  Dashed:  r-method 
using  two  relaxation  mechanisms.  Dotted:  Fade  approximant  method  using  two  relaxation 
mechanisms.  Dash-dotted:  Single  relaxation  mechanism  approximant. 
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t(s) 

FIG.  6.  Seismograms  collected  at  3000  m  offset  in  the  model  with  a  constant  Q  of 
approximation  of  20  between  2  and  25  Hz.  Solid:  Analytical  solution.  Dashed:  r-method 
using  two  relaxation  mechanisms.  Dotted:  Fade  approximant  method  using  two  relaxation 
mechanisms.  Dash-dotted:  Single  relaxation  mechanism  approximant. 
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FIG.  7.  Seismograms  collected  at  6000  m  offset  in  the  model  with  a  constant  Q  of 


approximation  of  20  between  2  and  25  Hz.  Solid:  Analytical  solution.  Dashed:  r-method 
using  two  relaxation  mechanisms.  Dotted:  Fade  approximant  method  using  two  relaxation 
mechanisms.  Dash-dotted:  Single  relaxation  mechanism  approximant. 
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FIG.  8.  Seismograms  collected  at  3000  m  offset  in  the  model  with  a  constant  Q  of 
approximation  of  20  between  2  and  25  Hz.  Solid:  Analytical  solution.  Dashed:  r-method 
using  five  relaxation  mechanisms.  Dotted:  Fade  approximant  method  using  five  relaxation 
mechanisms.  Dash-dotted:  Single  relaxation  mechanism  approximant. 
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FIG.  9.  Seismograms  collected  at  6000  m  offset  in  the  model  with  a  constant  Q  of 
approximation  of  20  between  2  and  25  Hz.  Solid:  Analytical  solution.  Dashed:  r-method 
using  five  relaxation  mechanisms.  Dotted:  Fade  approximant  method  using  five  relaxation 
mechanisms.  Dash-dotted:  Single  relaxation  mechanism  approximant. 
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